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ABSTRACT 


A  factor  of  safety  (FS)  method  for  quantitative  estimates  of  grid/time 
uncertainties  for  CFD  solutions  is  derived  to  remove  the  deficiencies  of  GCI,  corrected 
GCIc,  and  correction  factor  (CF)  methods,  i.e.,  unreasonably  small  uncertainty  when  CF 
>  1  (estimated  order  of  accuracy  greater  than  theoretical)  and  lack  of  statistical  analysis  to 
prove  95%  confidence  for  the  estimated  uncertainties  to  bound  the  true  error.  The 
approach  follows  the  CF  method  but  reflects  the  uncertainty  instead  of  FS  for  CF  <  1  for 
CF  >  1  (CF  =  1  is  asymptotic  range).  FS  at  CF  =  0  and  1  are  determined  by  reliability  and 
lower  band  of  the  confidence  interval  of  the  true  mean  based  on  statistical  analysis  using 
a  large  sample  of  analytical/numerical  benchmarks  covering  17  studies,  96  variables  and 
304  individual  grid  triplets.  Only  the  FS  method  provides  95%  confidence  that  the  actual 
factor  of  safety  FSa  >  1  for  the  304  grid  convergence  studies:  confidence  intervals  are 
86.2%,  92.1%,  91.5%,  and  95.7%  for  GCI,  GCIc,  CF,  and  FS.  For  20%  of  the  data  1.1  < 
CF  <  2.0,  GCI,  GCIc,  and  CF  methods  fail  as  only  47.4%,  71.2%,  72.9%  confidence 
intervals  are  achieved,  whereas  89.3%  is  achieved  for  the  FS  method.  Only  the  FS 
method  has  95%  confidence  the  lower  band  of  the  confidence  interval  for  FSa  is  larger 
than  1.2  for  different  studies,  variables,  ranges  of  CF,  and  single  CF  values  where 
multiple  FSa  are  available. 


IV 


I.  INTRODUCTION 


Current  quantitative  error  estimates  for  grid  size  and  time  convergence  are  based 
on  Richardson  extrapolation  (RE),  i.e.  the  error  is  expanded  in  a  power  series  expansion 
with  integer  powers  of  grid  spacing/time-step  as  a  finite  sum.  It  is  common  practice  to 
retain  only  the  1st  tenn  of  the  series  assuming  the  solutions  are  in  the  asymptotic  range 
(AR),  which  leads  to  a  grid  triplet  study.  The  grid  convergence  index  (GCI)  derived  by 
Roache  [1]  can  be  used  to  estimate  the  uncertainties  due  to  grid/time  errors  and  is  widely 
used  and  recommended  by  ASME  [2]  and  AIAA  [3].  Stem  et  al.  derived  the  correction 
factor  (CF)  method  [4]  with  improvements  made  by  Wilson  et  al.  [5].  They  introduced  a 
variable  factor  of  safety  (FS)  as  a  function  of  the  CF  so  that  FS  increases  linearly  with 
distance  of  solutions  from  the  AR.  CF  is  defined  by: 


,«Pre 


CF  = 


-1 


,,Pth 


(1) 


where  pR,:  is  the  estimated  order  of  accuracy  and  pth  is  the  theoretical  order  of  accuracy. 


CF  provides  an  important  metric  that  quantifies  the  distance  of  solutions  to  the  AR.  As 
shown  later,  CF  is  also  useful  for  statistical  analysis,  for  which  analytical  benchmark 
(AB)  or  numerical  benchmark  (NB)  data  can  be  organized  according  to  the  same  CF 
values.  The  CF  method  is  validated  for  CF  <  1  using  AB  including  one-dimensional  wave 
and  two-dimensional  Laplace  equations.  Since  the  AB  approach  AR  from  CF  <  1,  it  is 
assumed  that  the  FS  for  CF  >  1  is  obtained  by  reflecting  the  FS  for  CF  <  1  with  respect  to 
AR  (CF  =1)  assuming  FS  is  the  same  at  the  same  distance  from  the  AR. 

There  are  several  problems  of  using  RE.  As  shown  by  Stern  et  al.  [6],  it  is 
difficult  to  improve  the  accuracy  by  retaining  more  terms  in  the  power  series.  When 
solutions  are  not  in  the  AR,  multiple  grid  triplet  studies  often  show  non-smooth 
convergence,  i.e.,  pRE  approaches  pth  in  an  oscillatory  fashion  with  a  large  range  of 


values  [2].  RE  requires  at  least  3  systematic  high  quality  grids,  which  may  be  too 
expensive  for  industrial  applications.  The  grid  refinement  ratio  r  must  be  carefully 
selected,  r  cannot  be  too  large  as  the  grids  may  resolve  different  flow  physics.  Too  small 
values  of  r  (i.e.,  very  close  to  one)  are  also  undesirable  since  solution  changes  will  be 
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small  and  sensitivity  to  grid/time-step  may  be  difficult  to  identify  compared  to  iterative 
errors. 

The  non-smooth  grid  convergence  problem  may  be  resolved  using  the  least  square 
method  (LSM)  [7]  or  response  surface  method  (RSM)  [8],  which  requires  at  least  4 
solutions.  There  are  some  issues  of  implementing  the  LSM/RSM:  (1)  the  relationship 
between  the  LSM/RSM  and  the  individual  grid  triplet  studies  is  not  established;  (2)  it 
does  not  discriminate  between  converging  and  diverging  grid  studies  and  the  use  of 
diverging  solutions  is  not  well  founded;  (3)  the  requirement  of  at  least  four  gird  solutions 
is  often  too  expensive  for  industrial  applications;  (4)  all  solutions  are  required  to  be  in  the 
AR,  which  is  contradictory  to  the  use  of  non-smooth  and  non-monotonic  converged 
solutions;  and  (5)  it  introduces  additional  uncertainties  due  to  the  least-square  fit.  The 
expensive  cost  of  using  RE  may  be  resolved  using  the  single-grid  uncertainty  estimate 
approach  [9]  but  this  approach  is  not  well  developed  for  three-dimensional  problems  and 
does  not  provide  sensitivity  of  solutions  to  grid  size  and  time  step. 

GCI  and  CF  methods  have  two  deficiencies.  Uncertainty  estimates  when 
P re  >  Pth  (CF  >1)  are  unreasonably  small  in  comparison  to  those  with  the  same  distance 
to  AR  for  CF  <  1 .  This  is  due  to  the  fact  that  the  error  estimate  SRh:  for  the  former  is  much 

smaller  than  that  of  the  latter.  There  is  no  statistical  evidence  for  what  confidence  interval 
GCI  and  CF  can  actually  achieve.  It  is  claimed  by  Roache  [1]  and  the  ASME 
performance  test  codes  committee  PTC  61  [10]  that  a  95%  confidence  interval  is 
achieved  for  GCI  with  FS  =  1.25  based  on  over  500  demonstrated  cases  by  dozens  of 
groups.  However,  this  is  only  based  on  anecdotal  observations  as  no  statistical 
distributions  or  analysis  are  reported  in  the  references  cited  in  [1],  [10],  or  [1 1]  to  support 
this  claim. 

A  recent  study  by  Logan  and  Nitta  [8]  evaluated  10  different  methods  on 
uncertainty  estimates  using  limited  statistical  analysis  including  solution  verification 
method  reliability  Rsm  and  the  reduced  Chi-Square  A,2, .  A  higher  value  of  X2*  indicates 

that  a  solution  verification  method  is  too  conservative.  The  use  of  corrected  GCI  (GCIc) 
with  FS  =  1.25  is  much  closer  to  a  68%  confidence  estimate  than  95%.  There  is  no 
apparent  correlation  between  Rsm  and  X2*  ,  i.e.,  a  method  may  have  very  good  Rsm  but 
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Xvi  » 1 .  Since  this  study  only  has  3  structure  problems  with  18  individual  grid 

solutions,  it  was  suggested  that  a  larger  sample  set  with  the  number  of  grid  convergence 
sets  much  larger  than  100  is  needed  to  draw  any  general  conclusions. 

Two  other  recent  studies  [12,  13]  have  considered  the  use  of  different  uncertainty 
estimates  for  different  ranges  of  pRi:  for  monotonically  converged  solutions.  The  idea  is 

not  new  as  it  was  the  basis  for  the  CF  method  as  per  Stem  et  al.  [4],  although  the  CF 
method  was  not  referenced.  Eca  and  Hoekstra  [12]  extended  the  LSM  for  a  nominally 
second-order  accurate  method,  i.e.,  0  <  pm  <  0.95  (o<  CF  <o.3l),  0.95 <2.05 

( 0.3 1 <  CF  <  1.05),  and  pRE  >  2.05  {cf>  1.05).  According  to  the  authors,  the  estimates  were 

based  on  the  experience  obtained  in  a  variety  of  test  cases  and  the  suggestions  and 
comments  of  the  first  Workshop  on  CFD  Uncertainty  Analysis  [14].  Similar  ideas  were 
used  by  Rumsey  and  Thomas  [13]  who  provided  different  formulas  for  different  ranges 
of  p^  for  a  nominally  third-order  accurate  method,  i.e.,  0<  pm  <  0.95  (o<CF<o.l3), 
0.95  <  pRE  <  3.05  (0. 31  <  CF  <l.04),  and  pRE>  3.05  (cf>i.04).  Since  the  authors  did  not 
provide  detailed  derivation  or  validation  of  their  estimates  and  pm  ranges  it  is  difficult  to 

understand  their  differences.  Here  again  statistical  distributions  or  analysis  are  not 
reported. 

The  overall  objective  of  this  study  is  to  develop  a  FS  method  that  removes  the 
above  deficiencies  and  achieves  an  overall  95%  confidence  interval  for  the  estimated 
uncertainty  to  bound  the  true  error.  The  approach  is  to  extend  the  CF  method  by 
reflecting  uncertainty  not  FS  for  CF  <  1  for  CF  >  1  with  respect  to  the  AR.  The  FS 
method  provides  the  flexibility  of  specifying  FS  at  CF  =  0  and  CF  =  1,  i.e.,  FSo  and  FSi, 
which  will  be  determined  using  rigorous  verification  studies  covering  a  wide  range  of 
problems  that  have  either  AB  or  NB  data.  Statistical  analysis  includes  the  reliability  and 
confidence  interval  for  the  mean  of  the  actual  factor  of  safety  (FSa).  The  FS  method  will 
be  compared  with  GCI,  GCIc,  and  CF  methods.  It  should  be  noted  that  the  idea  of 
reflecting  uncertainty  not  FS  for  CF  <  1  for  CF  >  1  with  respect  to  the  AR  is  not  new  as  it 
was  first  reported  in  an  earlier  study  by  the  authors  [15],  However,  the  original  idea  was 
developed  only  logically  based  on  observations  that  uncertainty  estimates  for  CF  >  1  are 
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unreasonably  small  and  simply  tested  on  one  example  of  ship  hydrodynamics  application. 
It  was  not  validated  by  AB  and  NB  as  conducted  herein.  Additionally,  the  original  idea  is 
altered  so  that  users  can  flexibly  specify  FSo  and  FSi.  Due  to  these  facts,  the  present  work 
supersedes  the  previous  report  [15]. 

II.  ERROR  AND  UNCERTAINTY  ESTIMATE  USING  RICHARDSON 

EXTRAPOLATION 

Grid/time  convergence  studies  are  conducted  with  multiple  solutions  using 
systematically  refined  grid  sizes  or  time  steps.  First  r  for  grid/time  is  selected.  As  an 
example,  if  3,  2,  and  1  represent  the  coarse,  medium,  and  fine  grids  with  grid  spacing 
A*3,  Ax 2,  and  Axi,  respectively,  then 

Ax,  Ax, 

r  =  — 1  1  (2) 

Axj  Ax2 

Constant  r  is  not  required  [1]  but  simplifies  the  analysis  and  thus  used  herein.  Solution 
changes  s  for  medium-fine  and  coarse-medium  solutions  and  the  convergence  ratio  R  are 
defined  by 

£2\  ~  S2  ~  Sl 

£32=S,-S2  (3) 

R  —  ^21/^32 

When  0  <  R  <  1  solution  monotonic  convergence  is  achieved  and  generalized  RE  is  used 
to  estimate  order-of-accuracy  pRE ,  the  error  SjtF  and  numerical  benchmark  Sc-  The  error 

is  expanded  in  a  power  series  expansion  with  integer  powers  of  grid  spacing/time-step  as 
a  finite  sum.  The  accuracy  of  the  estimates  depends  on  how  many  terms  are  retained  in 
the  expansion,  the  magnitude  (importance)  of  the  higher-order  terms,  and  the  validity  of 
the  assumptions  made  in  RE  theory.  In  general,  2n+ 1  solutions  are  required  to  estimate 
Sc  and  the  first  n  terms  in  the  power  expansion.  When  solutions  are  close  to  the  AR,  grid 
spacing  is  sufficient  small  and  the  lower-order  terms  in  the  expansion  dominate.  If  only 
the  first  term  is  retained  (n=  1),  three  grid  solutions  are  needed  to  provide  estimates  for 

P RE  ’  ^RE  ’  Sc 
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^  (^32/^21) 

Pre~  ln(r) 

(4) 

8  =  £lx 

M  rP-  _! 

(5) 

Sc  =  s-sa 

(6) 

Eqn.  (5)  provides  the  error  estimate  with  both  sign  and  magnitude  but  not  an  “error  band” 
(uncertainty)  in  which  the  users  can  have  some  practical  level  of  confidence. 

Roache  proposed  the  GCI  method.  The  idea  behind  GCI  is  to  approximately  relate 
the  821  in  Eqn.  (5)  obtained  by  whatever  grid  convergence  study  (whatever  pRE  and  r)  to 
the  8  that  would  be  expected  from  a  grid  convergence  study  of  the  same  problem  with  the 
same  fine  grid  using  pRE  =  2  and  r  =  2,  i.e.  a  grid  doubling  with  a  2nd-order  method.  The 
relation  is  based  on  equality  of  the  error  estimates.  Given  an  821  from  an  actual  grid 
convergence  test,  the  GCI  is  derived  by  calculating  the  error  estimate  SRE  from  Eqns.  (5) 
and  (4),  then  calculating  an  equivalent  821  that  would  produce  approximately  the  same 
SKF  with  pRE  =  2  and  r  =  2.  The  absolute  value  of  that  equivalent  8  is  the  GCI  for  the  fine 
grid  solution,  which  is  conveniently  expressed  as  [1] 

t/^=ra^ri=ra|‘^1  <7) 


FS  can  thus  be  regarded  as  a  “factor  of  safety”  over  the  RE  error  estimate  8m .  Roache 
suggested  FS  =  1.25  for  systematic  grid-triplet  studies  using  RE  estimate  for  order  of 
accuracy  pRE  and  FS  =  3  for  2-grid  sensitivity  studies  using  theoretical  estimate  for  order 

of  accuracy  pth .  In  Fig.  1,  FS  =  1.25  is  plotted  versus  CF,  which  suggests  constant  FS  for 
all  ranges  of  CF  values. 

Fogan  and  Nitta  [8]  also  used  GCIc,  which  modifies  Eqn.  (7)  by  discarding  p m 
and  use  pth  when  pRK  >  pth  .  As  can  be  easily  shown,  GCIc  is  equivalent  to  multiply  FS 
by  CF  when  CF  >  1 ,  which  leads  to 


U, 


GCI 


FSxCF\Srf 


CF  <  1 
CF  >  1 


(8) 
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The  CF  method  was  derived  based  on  AB  for  CF  <  1,  which  shows  that  an 


improved  error  estimate  is  to  multiply  Eqn.  (5)  by  CF 


CF5m  =  CF 


"21 


„Pre 


(9) 


The  use  of  Eqn.  (9)  will  replace  p ^  using  pth .  However,  pm  is  not  discarded  but 


included  in  CF.  The  uncertainty  is  estimated  by  the  sum  of  the  absolute  value  of  the 
corrected  error  estimate  from  RE  and  the  absolute  value  of  the  amount  of  the  correction 

U  =  FS(CF)\SJ  =  I  CF8re  I  + 1(1  -  CF)  Sm  I  (10) 


As  criticized  by  Roache  [16],  Eqn.  (10)  is  deficient  for  CF  <  1  in  only  providing  50% 
confidence  level.  Wilson  et  al.  [5]  revised  Eqn.  (10)  for  proper  behavior  for  CF  <  1  by 
increasing  FS  for  decreasing  CF,  i.e.  FS  is  a  function  of  CF.  Then  the  variable  FS  for  CF 
<  1  is  reflected  for  CF  >  1  with  respect  to  the  AR  (CF  =  1).  Additionally,  they  modified 
the  formula  for  proper  behavior  for  CF  =  1  by  providing  10%  FS  in  the  limit  CF  =  1 
while  smoothly  merging  with  previous  CF  method  uncertainty  estimates  for 
1 1  -  CF  |>  0. 125  for  the  uncorrected  solutions. 


U  =  FS{CF)Sre 


[9.6(1-CF)2+1.1 
[2|l-CF|  +  l]|^| 


0.875  <CF<  1.125 


0<CF<  0.875  or  CF  >  1.125 


(ID 


As  shown  by  Wilson  et  al.  [5,  17]  and  Fig.  1,  the  CF  method  is  equivalent  to  the  GCI,  but 
with  a  variable  FS.  The  variable  FS  has  the  “common-sense”  advantage  in  providing  a 
quantitative  metric  to  determine  proximity  of  the  solutions  to  the  AR  and  approximately 
accounts  for  the  effects  of  higher-order  RE  terms.  The  CF  method  has  been  used  in  ship 
hydrodynamics  CFD  workshops. 


III.  FACTOR  OF  SAFETY  METHOD 

The  overall  approach  of  the  FS  method  follows  the  CF  method  [4,  5]  with  a  linearly 
varying  FS  for  CF  <  1 .  However,  it  has  two  improvements  on  top  of  the  CF  method  by 
introducing  two  flexible  parameters,  FSo  and  FSi.  Additionally,  it  reflects  the  uncertainty 
rather  than  the  FS  itself  for  CF  <1  for  CF  >  1  with  respect  to  the  distance  from  CF  =  1. 
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Recommended  values  of  FSo  and  FSi  are  determined  using  large  sample  of  AB  or 
NB  and  statistical  analysis  compared  to  only  AB  and  no  statistical  analysis  for  the  CF 
method.  The  procedure  is  to  minimize  the  two  parameters  until  two  criteria  are  met:  (1) 
overall  at  least  95%  confidence  interval  is  achieved  for  the  FSa  is  larger  than  1;  (2)  at 
least  95%  confidence  that  the  lower  band  of  the  confidence  interval  for  the  true  mean  FSa 
is  larger  than  1.2. 

Following  Stem  et  al.  [4]  with  improvements  by  Wilson  et  al.  [5],  the  FS  method 
first  assumes  that  uncertainty  estimate  has  the  same  form  as  Eqns.  (10)  and  (11)  for 

CF  <  1 ,  i.e.  U  =  FS  (CF)  <5^1  and  FS  is  a  linear  function  of  CF  that  increases  with  the 

distance  from  CF  =  1.  However,  FSi  and  FSo  are  multiplied  before  the  corrected  estimate 
from  RE  and  the  absolute  value  of  the  amount  of  the  correction,  respectively: 


U  =  FS,\  CFS^,  |  +FS„  \(\-CF)Sre\ 

(12) 

(7  =  [FS,  |  CF  |  +FS,  1  -CF  |] |  <5^.  | 

(13) 

~FS(CF) 


For  0  <  CF  <  1 ,  the  above  formula  becomes 

U  =  [FS{CF  +  FS0(l  -  CF)]  \SRE  |  0  <  CF  <  1  (14) 

To  overcome  the  too  small  uncertainty  estimate  for  CF  >  1 ,  the  uncertainty  instead  of  the 
FS  in  CF  <  1  is  reflected  for  CF  >  1  with  respect  to  the  distance  from  the  AR  [15].  First, 

r Pm  in  Eqn.  (5)  is  re-expressed  based  on  the  definition  of  CF: 

rpRE  =CF(rp*-\)  +  \  (15) 

Second,  Eqn.  (5)  is  substituted  into  Eqn.  (14)  with  the  use  of  Eqn.  (15),  which  results  in 
an  alternative  form  of  U 


U  =  [FS]CF  +  FS0(\-CF)\ 


0  <  CF  <  1 


(16) 


CF  in  the  above  equation  is  replaced  by  2 -CF .  Thus  for  the  same  r,  pth,  and  e2{ , 
Equation  (16)  becomes: 
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CF 


1  <  CF  <  2 


(17) 


[FSl(2-CF)  +  FS0(CF-\)]\SRE 


u  = - 

2  -CF 

In  summary,  uncertainty  for  uncorrected  solutions  using  the  FS  method  is: 

[FSXCF  +  FS0(\  -  CF)]  \8rf  |  0  <  CF  <  1 


U  =  \ 


CF 


2 -CF 


[FSl(2-CF)  +  FS0(CF-\)]SRE 


1  <  CF  <  2 


(18) 


Compared  to  the  CF  method,  the  FS  method  introduces  an  additional  tenn  CF j{  2  -  CF) 

to  compute  U  for  CF  >  1.  When  CF  increases  from  1  to  2,  this  tenn  increases  rapidly 
from  1  to  infinity,  which  amplifies  FS  when  solutions  are  further  away  from  the  AR.  The 
FS  method  is  only  applicable  for  0  <  CF  <  2.  CF  =  0  is  the  border  of  convergence  and 
divergence  such  that  grid  errors/uncertainties  are  infinite  due  to  infinite  8m  as  a  result  of 


pRF  =0,  i.e.  solution  changes  for  the  medium  and  fine  grids  are  equal  to  those  for  the 


coarse  and  medium  grids.  For  CF  >  2,  solutions  are  too  far  from  the  AR  and  also 
regarded  as  divergent.  Figure  1  compares  the  theoretical  FS  predicted  by  GCI,  GCIc,  CF, 
and  FS  methods.  FS  =  1.0  suggests  that  the  design  meets  but  does  not  exceed  the 
minimum  requirements  with  no  room  for  variation  nor  error  while  too  high  FS  results  in 
excessive  weight  and/or  cost. 

The  FS  method  provides  users  flexibility  in  specifying  FSo  and  FSi.  To  determine 
the  “optimal”  values  for  FSo  and  FSi  and  validate  the  FS  method  requires  “iterations” 
based  on  statistical  analysis  of  a  large  sample  of  data  (see  Section  4).  The  iterations 
outcome  are  FSo=  2  and  FSi  =  1.25: 


U  = 


(2-0.75CF)  |  SRE  | 
CF(0.5  +  0.75CF) 


0  <  CF  <  1 


'RE 


(19) 


2 -CF 


1  <  CF  <  2 


Figure  1.  Factors  of  safety  for  GCI,  GCIc,  CF,  and  FS  verification  methods. 


IV.  STATISTICAL  ANALYSIS 


Reliability  R1  is  defined  as 

_  number  of  studies  FS A  >  1 


(20) 


Total  number  of  studies 
where  FSa  is  defined  as  the  ratio  of  the  uncertainty  estimate  and  the  true  error  Ei  between 
the  fine  grid  solution  and  the  AB/NB  solution  Sc : 


fsa=T 

E, 


Ey  = 


S\-Sc 


Sr 


(21) 

(22) 


R™=1 -(Estimated  fraction  for  FSa>1  -  Expected  fraction  for  FSa  >11  used  in  [8]  is 
deficient  as  it  does  not  discriminate  estimated  fractions  with  the  same  distance  below  or 


1  R  hereinafter  not  to  be  confused  with  R  used  before  for  convergence  ratio 
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above  the  expected  fraction  (e.g.,  91%  estimated  fraction  inside  has  the  same  Rsm  as  99% 
estimated  fraction  inside). 

The  confidence  of  mean  analysis  is  based  on  the  methodology  and  procedures 
summarized  in  [18].  If  X.  (i  =  1  ,N)  is  any  individual  sample  population/distribution,  the 


sample  mean  of  X  is 


_  1  N 

Ntt  1 

where  N  is  the  sample  size.  The  standard  deviation  of  the  sample  is  defined  by 


(23) 


5 


x , 


N  n 

i= 1 _ 

N- 1 


(24) 


and  is  related  to  the  standard  deviation  for  the  mean  S—  by 


sx 

Sj'~7n 


(25) 


To  account  for  the  effect  of  small  number  of  the  sampling  data,  student  t- 
distribution  is  applied  where  the  confidence  interval  for  the  mean  is  defined  as: 

tSx 

(26) 


k=tS-*  =  4N 


The  confidence  interval  of  mean  is  that  at  least  95%  confidence  that  true  mean  f.1 


is  bounded  by  X-k  and  X  +  k 

p(X-k<ju<X  +  k)>  0.95  (27) 

The  “outliers”  for  all  the  statistical  analysis  are  identified  using  the  Peirce’s 
Criterion  as  summarized  in  [19].  Compared  to  the  Chauvenet’s  criterion,  Peirce’s 
criterion  is  more  rigorous,  does  not  make  an  arbitrary  assumption  concerning  the 
rejection  of  data,  and  theoretically  accounts  for  the  case  where  there  is  more  than  one 
suspect  data. 

The  sample  populations  are  constructed  based  on  statistical  analysis  of  AB  and 
NB  selected  from  17  different  studies  covering  9  fluids,  5  thermal  and  3  structure 
problems  with  additional  information  summarized  in  Table  1.  The  17  studies  were  solved 
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Table  1.  Verification  studies 


# 

Case 

Conditions 

Verification  variables 

Grids 

Outlier 

AB/NB 

1 

1D£  wave  [6] 

- 

Wave  profile 

10  grids 

0 

AB 

2 

2D  Laplace  [5,7] 

- 

Arbitrary  function 

6  grids 

0 

AB 

3 

2D  driven  cavity 
[20] 

Re=1000 

Maximum/minimum 

Of  stream-function, 
vorticity 

4  grids 

1 

NB  [21] 

4 

2D  natural 
convection  flows  in 
square  cavities  [22] 

Ra=104 

Max.  or  monitored 
Velocity,  location, 
temperature,  and  Nu 

5  grids 

2 

NB 

5 

2D  natural 
convection  flows  in 
square  cavities  [22] 

Ra=105 

Max.  or  monitored 
Velocity,  location, 
temperature,  and  Nu 

6  grids 

1 

NB 

6 

2D  natural 
convection  flows  in 
square  cavities  [22] 

Ra=106 

Max.  or  monitored 
Velocity,  location, 
temperature,  and  Nu 

7  grids 

0 

NB 

7 

Backward-facing 
step  [23] 

Re=l.  5x10s 

Reattachment  length, 
velocity 

7  grids 

1 

NB 

T24,  25] 

8 

2D  driven  cavity 
126,  27] 

Re=100 

Velocity 

5  grids 

0 

NB 

9 

2D  driven  cavity 
126,  27] 

Re=1000 

Velocity 

5  grids 

0 

NB 

10 

3D  cubic  cavity 
[26,27] 

Re=100 

Velocity 

4  grids 

3 

NB 

11 

Axisymmetric 
turbulent  flow 
through  a  valve 
[28,  29] 

Re=105 

Velocity,  TKE,  e 

5  grids 

4 

NB 

12 

ID  steady-state 
convection- 
diffusion  [11] 

Pe=l  and 
Pe=10 

Arbitrary  function 

6  grids 

0 

AB 

13 

Isothermal  cylinder 
enclosed  by  a 
square  duct 
[27,  28,  29] 

Ra=106, 

Pi- 10 

Velocity,  temperature 

5  grids 

1 

NB 

14 

Premixed 

methane/air  laminar 
flat  flame  on  a 
perforated  burner 
[28,30,311 

Inlet 

temperature 

298.2K 

Velocities,  temperature 

7  grids 

2 

NB 

15 

Data  for  “exact” 
grid  convergence 
set  [8] 

contrived 

- 

7  grids 

0 

AB 

16 

Bean  bending 
problem  for  2nd 
series  [8] 

2nd  series 

Beam  bending  stress, 
Beam  end  deflection 

7  grids  with  5 
systematically 
refined 

3 

AB 

17 

Beam  bending 
problem  for  3rd 
series  [8] 

3rd  series 

Beam  bending  stress, 
Beam  end  deflection 

4  grids 

0 

AB 

1 D  stands  for  “dimensional” 


on  four,  five,  or  six  refinement  levels  using  a  systematic  grid  refinement  ratio  2  (r  =  2) 
for  different  kinds  of  flows,  using  different  grids  and  different  numerical  schemes.  The 
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17  studies  constitute  96  variables  (velocity,  locations  max/min,  temperature,  turbulent 
kinetic  energy,  stream  function,  beam  bending  stress,  and  beam  end  deflection)  of  which 
either  one  or  multiple  grid  convergence  studies  are  available,  which  gives  a  total  of  304 
individual  monotonically  converged  grid  studies.  304  grid  studies  have  304  CF  values  of 
which  at  least  5  FSa  are  available  at  16  different  CF.  Details  about  the  numerical 
schemes,  solver,  and  boundary  conditions  can  be  found  in  the  corresponding  references. 

The  statistical  analysis  is  performed  in  different  ways  resulting  in  a  total  of  24 

different  distributions.  For  each  distribution,  variables  R,  X,  S-,  and  X-tS-  are 

calculated.  As  summarized  in  Table  2,  the  1st  and  2nd  distributions  are  for  the  17  means  of 
studies  and  the  means  of  the  96  variables  of  the  17  studies,  respectively.  Then  statistical 
analysis  is  conducted  for  different  ranges  of  CF  where  an  averaging  process  is  performed 
using  a  tolerance  ACF  =  0.01,  i.e.,  CF  values  are  regarded  to  be  the  same  if  the  difference 
between  two  CF  values  is  less  than  0.01.  A  smaller  value  with  ACF  =  0.001  is  evaluated 
and  no  significant  effects  on  the  statistics  are  observed.  After  the  averaging  process,  the 
3ld  to  8th  distributions  are  analyzed  for  different  ranges  of  CF  values,  i.e.,  (0-2),  (0-0.4), 
(0.4-0. 9),  (0. 9-1.1),  (1.1-1. 6),  and  (1. 6-2.0),  for  which  the  results  the  summarized  in 
Table  3.  Table  4  shows  the  9th  to  24th  distributions  that  are  at  the  16  CF  values  where 
there  are  more  than  4  monotonic  converged  grid  studies  and  thus  more  than  4  FSa  for 
each  CF.  The  numbers  of  outliers  for  each  study  and  the  16  CF  where  multiple  FSa  (>4) 
are  available  are  summarized  in  Table  1  and  Table  4,  respectively. 


Table  2.  Mean,  standard  deviation  of  the  mean,  and  reliability  excluding  outliers  based  on 
different  studies 


Factor  of  safety 

#  points 

Statistics 

GCI 

GCIc 

CF 

FS 

t 

Means  of 
studies 

17 

X 

1.61 

1.72 

2.27 

3.25 

1.75 

% 

0.16 

0.16 

0.31 

0.81 

X-tSJ 

1.33 

1.44 

1.73 

1.83 

R(%) 

82.35 

94.12 

94.12 

94.12 

Means  of  all 
variables  of  1 7 
studies 

96 

X 

1.63 

1.73 

2.24 

2.74 

1.67 

0.09 

0.09 

0.19 

0.28 

X-tST 

1.48 

1.58 

1.92 

2.27 

R  (%) 

88.54 

91.67 

91.67 

95.83 
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Table  3.  Mean,  standard  deviation  of  the  mean,  and  reliability  excluding  outliers  based  on 


different  ranges  of  correction  factors  using  non-averaged  factor  of  safety 


CF 

#  points 

Statistics 

GCI 

O 

O 

1— 1 
n 

CF 

FS 

t 

0-2 

304 

X 

1.69 

1.76 

2.34 

2.60 

1.65 

% 

0.11 

0.11 

0.27 

0.24 

1.51 

1.58 

1.89 

2.20 

R  (%) 

86.18 

92.11 

91.45 

95.72 

0-0.4 

27 

(8.88%) 

X 

4.77 

4.77 

10.20 

7.16 

1.71 

% 

1.06 

1.06 

2.54 

1.70 

X-tSj 

2.96 

2.96 

5.86 

4.25 

R  (%) 

96.30 

96.30 

100 

96.30 

0.4-0. 9 

96 

(31.58%) 

X 

1.74 

1.74 

2.24 

2.06 

1.67 

% 

0.08 

0.08 

0.14 

0.11 

1.61 

1.61 

2.01 

1.88 

R  (%) 

93.75 

93.75 

95.83 

95.83 

0.9-1. 1 

125 

(41.11%) 

X 

1.29 

1.30 

1.15 

1.34 

1.98 

0.02 

0.02 

0.02 

0.02 

1.25 

1.26 

1.11 

1.30 

R  (%) 

96.0 

96.0 

94.4 

98.4 

1. 1-1.6 

48 

(15.79%) 

X 

1.02 

1.29 

1.24 

2.16 

1.68 

% 

0.06 

0.07 

0.07 

0.15 

x-ts T 

0.92 

1.17 

1.12 

1.91 

R  (%) 

47.92 

77.08 

70.83 

87.5 

1. 6-2.0 

8 

(2.63%) 

X 

0.97 

1.74 

2.01 

16.09 

1.90 

% 

0.14 

0.25 

0.29 

3.80 

x~tsJ 

0.70 

1.27 

1.46 

8.87 

R  (%) 

37.5 

37.5 

87.5 

100 

Table  4.  Statistics  excluding  outliers  at  16  different  correction  factors 


CF 

#  points 

#outliers 

Statistics 

GCI 

GCIc 

CF 

FS 

t 

0.625 

4 

1 

X 

1.69 

1.69 

2.36 

2.07 

2.35 

% 

0.09 

0.09 

0.12 

0.11 

x~tSJ 

1.48 

1.48 

2.07 

1.81 

R  (%) 

100 

100 

100 

100 

0.675 

9 

0 

X 

1.90 

1.90 

2.51 

2.27 

1.86 

0.12 

0.12 

0.15 

0.14 

X~*x 

1.69 

1.69 

2.23 

2.02 

R  (%) 

100 

100 

100 

100 
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X 

1.26 

1.27 

1.11 

1.27 

1.005 

50 

5 

% 

0.025 

0.025 

0.022 

0.025 

1.68 

X-tSj 

1.22 

1.22 

1.08 

1.23 

R(%) 

96 

96 

94 

96 

X 

1.24 

1.26 

1.09 

1.29 

1.015 

6 

0 

% 

0.008 

0.008 

0.007 

0.010 

2.02 

X-tS^ 

1.23 

1.24 

1.08 

1.27 

R(%) 

100 

100 

100 

100 

X 

1.72 

1.89 

1.64 

2.21 

1.095 

5 

0 

h 

0.17 

0.19 

0.16 

0.22 

2.13 

X-tS* 

1.36 

1.49 

1.30 

1.75 

R(%) 

100 

100 

100 

100 

X 

1.19 

1.36 

1.23 

1.73 

1.145 

13 

1 

% 

0.11 

0.12 

0.11 

0.15 

1.78 

X-tS^ 

1.00 

1.15 

1.03 

1.45 

R(%) 

69.23 

76.92 

69.23 

92.31 

X 

1.09 

1.42 

1.39 

2.37 

1.295 

8 

2 

0.05 

0.06 

0.06 

0.11 

1.90 

X-tSx 

1.00 

1.29 

1.27 

2.17 

R  (%) 

75 

100 

100 

100 

Determination  of  the  optimal  values  for  FSo  and  FSi  require  iterations  since 
computation  of  the  mean,  standard  deviation,  and  k  require  the  values  of  FSo  and  FSi 
known  as  a  priori.  Two  criteria  are  used  to  determine  if  the  selected  values  are 
acceptable.  The  first  criterion  is  that  for  the  total  304  grid  convergence  studies  at  least 
95%  confidence  is  achieved  for  FSa>  1,  i.e., 

R  >  95%  (28) 

R  is  also  used  to  estimate  the  performance  of  different  verification  methods  for  different 
studies,  different  regions  of  CF  values,  and  at  the  16  different  CF  values  with  multiple 
solutions.  The  2nd  criterion  is  that  there  will  be  at  least  95%  confidence  that  the  lower 
band  of  the  confidence  interval  for  the  mean  FSa  is  larger  than  FSA,min-  The  appropriate 
values  of  FSA,min  are  often  determined  based  on  risks/reliability,  accuracy,  and  cost,  etc. 
and  show  a  large  variations  from  1.2  to  5  for  different  practical  applications.  For 
example,  1.2  is  used  for  new  bridges  and  road  marks,  3.0  is  used  for  automobiles,  and  3.5 
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to  4.0  are  used  for  pressure  vessels,  etc.  FSa, min  is  thus  chosen  to  be  1.2  in  this  study 


p(l.2  <  X—k  <  ju<  X +k^>  0.95 

(29) 

It  is  obvious  that  when  Eqn.  (29)  is  satisfied 

p{x>  1.2)  >0.95 

(30) 

P(ju>  1.2)  >0.95 

(31) 

When  the  two  criteria  are  met,  a  minimum  set  of  values  are  accepted.  In  this  study,  (FSo  = 
3,  FSi  =  2),  (FS0  =  3,  FSi  =  1.5),  (FS0=  2,  FSi  =  1.5)  and  (FS0=  2,  FSi  =  1.25)  are  found 
to  meet  the  two  criteria.  Since  the  first  three  sets  of  data  are  too  conservative  (i.e,  provide 
too  large  FSa),  (FSo=  2,  FSi  =  1.25)  is  determined  to  be  the  recommended  values,  which 
are  used  and  compared  with  GCI,  GCIc,  and  CF  methods  for  the  same  sample  of  data  to 
evaluate  the  mean,  standard  deviation  for  the  mean,  lower  band  value  for  the  confidence 
interval  of  the  mean  FSa,  and  the  reliability  based  on  different  studies  and  different 
ranges  of  CF  values. 

It  should  be  noted  that  statistical  analysis  is  perfonned  based  on  the  following 
fact.  The  error  and  uncertainty  estimates  are  systematic,  but  since  the  real  error  Ei  is 
random,  the  ratio  between  them  (i.e.,  FSa)  is  randomly  distributed.  In  other  words,  it  is 
assumed  that  there  are  no  correlated  biased  errors  between  different  studies.  Since  FSa  is 

randomly  distributed,  the  confidence  interval  for  the  mean  reveals  how  close  X  is  to  |l 
for  FSa- 


V,  ANALYTICAL  AND  NUMERICAL  BENCHMARK  DATA 

All  the  17  studies  summarized  in  Table  1  used  the  same  grid  refinement  ratio  2  to 
generate  the  grids  systematically.  For  each  study,  solutions  on  multiple  grids  are 
extracted  and  convergence  studies  are  performed  for  each  variable.  Only  monotonically 
converged  grid  studies  with  0  <  CF  <  2  are  used,  which  constitute  about  90%  of  all  the 
grid  convergence  studies.  The  other  10%  with  oscillatory  convergence,  monotonic 
divergence,  and  oscillatory  divergence  are  discarded,  pre,  Sre,  and  Sc  are  evaluated  using 
Eqns.  (4),  (5),  and  (6),  respectively. 
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For  each  study,  two  types  of  figures  are  plotted  and  analyzed.  To  investigate  the 
convergence  characteristics  of  each  variable,  |U|,  actual  error  |Ei|,  and  CF  are  graphed 
versus  the  ratio  of  grid  spacing  Ax/Axfinest.  Additionally,  FSa  is  plotted  versus  CF  to  show 
the  actual  behavior  of  each  method  for  different  ranges  of  CF  values,  which  enables  the 
evaluation  of  different  methods  on  whether  the  uncertainty  estimate  bound  the  true  error 
(|U|  >|Ei|).  Details  of  these  figures  and  discussions  are  presented  in  Appendix  1. 

37  of  the  96  variables  only  have  two  grid  triplet  studies,  which  prohibit  studies  of 
the  convergence  characteristics  of  those  variables.  However,  the  errors  do  show 
decreasing  values  when  the  grids  are  refined.  Of  the  remaining  59  variables  that  have 
multiple  grid  triplet  studies,  ideal  convergence  is  only  achieved  for  6  of  the  variables 
such  as  the  one-dimensional  wave  equation,  i.e.,  errors  predicted  by  the  four  methods 
monotonically  decrease  continuously  as  the  solutions  are  approaching  the  AR  when  grids 
are  refined.  The  other  variables  oscillatorially  approach  the  AR  when  the  grids  are 
refined.  FSa  of  all  variables  for  each  study  are  also  examined  for  different  CF  values, 
which  will  be  discussed  in  more  details  in  the  next  section. 

VI.  STATISTICAL  ANALYSIS  RESULTS 

Figure  2  shows  the  FSa  without  averaging  and  without  excluding  outliers  using 
the  four  different  methods  for  the  total  304  verification  studies.  There  are  some  cases  at 
CF  around  0.5,  0.65,  0.8,  1,  and  1.4  that  none  of  the  methods  achieves  FSa>  1-  As  shown 
later,  some  of  these  points  are  outliers  but  the  rest  of  them  are  not.  Compared  to  the  CF 
method,  GCI  predicts  a  lower  and  higher  FSa  near  CF  =  0  and  CF  =  1,  respectively, 
which  is  consistent  with  the  theoretical  FS  distribution  shown  in  Figure  1.  Compared  to 
GCI,  GCIc  predicts  a  little  higher  FSa  for  CF  >  1  but  still  not  high  enough  to  bound  most 
of  the  errors.  Both  GCI  and  the  CF  methods  do  not  provide  sufficient  large  FSa  for  CF  > 
1,  which  is  corrected  by  the  FS  method  as  it  shows  a  nearly  “symmetric”  distribution  of 
FSa  respect  to  CF  =  1 . 

Figure  3  shows  FSa  with/without  averaging  using  ACF  =  0.01  and  X±tS-  for 

CF  where  at  least  four  FSa  are  available  excluding  outliers  using  the  four  methods.  At  the 
CF  where  multiple  FSa  are  available,  the  number  of  sampled  data  from  left  to  right:  4,  9, 
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7,  14,  12,  4,  9,  7,  4,  7,  7,  50,  6,  5,  13,  8  points.  It  is  observed  that  the  averaging  itself  only 
helps  to  remove  a  few  “unbounded”  points  but  still  cannot  drag  all  the  FSa  to  be  larger 
than  1  for  all  the  methods.  For  X-tS-,  GCI  has  FSa  <  1  for  most  CF  >  1.  GCIc 

increases  the  FSa  magnitude  for  CF  >  1  but  still  not  large  enough  to  be  greater  than  1 .2 
for  all  CF  where  multiple  FSa  are  available.  Compared  to  GCI  and  GCIc,  CF  predicts 
lower  FS  near  CF  =  1  and  similar  magnitude  of  FS  for  CF  >  1 .  Only  the  FS  method 
bounds  the  largest  fraction  of  the  error  and  provides  a  minimum  FSa  >1.2  for  all  CF 
where  multiple  FSa  are  available. 


(c) 

Figure  2  Actual  factor  of  safety  without  averaging  and  without  excluding  outliers: 
(a)  GCI  and  GCIc,  (b)  CF  method,  (c)  FS  method. 
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Figure  3  Actual  factor  of  safety  with/without  averaging  using  ACF  =  0.01  and  X  ±  tS  x 

for  CF  where  multiple  FSa  are  available  (number  of  sampled  data  from  left  to 
right:  4,  9,  7,  14,  12,  4,  9,  7,  4,  7,  7,  50,  6,  5,  13,  8  points,  respectively) 
excluding  outliers:  (a)  GCI,  (b)  GCIc,  (c)  CF  method,  (d)  FS  method. 

Figure  4  shows  the  standard  deviation  based  on  the  mean  (Sx  jx)  for  the 
different  methods  excluding  outliers  at  the  16  CF  where  multiple  FSA  are  available. 
Sx  j X  is  the  same  for  different  methods  as  dividing  by  the  mean  cancels  the  different 

factors  of  safety  used  in  the  formula.  Overall  Sx  j X  shows  strong  oscillations  with 
magnitude  decreasing  when  approaching  the  AR.  This  is  consistent  with  the  facts  that 
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larger  number  of  verification  studies  are  available  close  to  the  AR  than  those  far  away 
from  AR,  which  will  lead  to  smaller  Sx  . 


Figure  4  Standard  deviation  based  on  mean  for  different  methods  excluding  outliers  at  the 
16  CF  where  multiple  FSa  are  available. 

Tables  2  and  3  show  the  statistics  for  the  four  methods,  either  based  on  different 
studies  and  variables  (Table  2)  or  based  on  different  ranges  of  CF  values  (Table  3).  The 
outliers  have  been  removed.  Based  on  17  studies,  none  of  the  methods  achieves  95%  (R  > 
95%)  likely  due  to  the  small  population  number.  GCI  has  the  smallest  R  =  82.4%  while 

the  other  three  methods  have  higher  R  =  94.1%.  X  of  FSa  for  the  17  studies  are 
increasing  from  1.61  for  GCI  to  3.25  for  FS.  S-  is  increasing  from  0.16  for  GCI  to  0.81 

for  FS  method.  Based  on  all  96  variables  considered,  only  the  FS  method  has  R  >  95% 
while  R  for  other  methods  are:  R  (GCI)  <  R  (GCIc)  =  R  (CF)  <  R  (FS).  Compared  to 
values  using  17  studies,  GCI,  GCIc,  and  CF  methods  show  no  significant  changes  of  X 
but  significant  lower  S-.  A  and  S-  of  the  FS  method  decrease  to  2.74  and  0.28, 
respectively.  For  all  four  methods,  X-tSx  >1.2  is  met  based  on  either  the  17  studies  or  96 
variables. 

FSa  presented  in  Table  2  based  on  different  studies  and  variables  show  larger 
mean  FSa  of  FS  method  than  any  other  method.  This  may  be  misleading  as  it  seems  that 
the  higher  reliability  using  the  FS  method  is  due  to  the  larger  FS.  In  fact,  statistical 
analysis  using  different  studies  involves  averaging  of  different  verification  studies  that 
have  0  <  CF  <  2 .  This  will  increase  the  overall  mean  significantly  as  the  FS  far  away 
from  the  AR  is  much  larger  than  FS  close  to  the  AR  for  the  CF  and  FS  methods.  To 
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illustrate  the  distribution  of  FSa  and  evaluate  the  effectiveness  of  different  methods  in 
different  ranges  of  CF,  Table  3  shows  statistics  for  different  ranges  of  CF  using  non- 
averaged  FSa  excluding  outliers.  For  the  whole  range  of  CF  (0-2):  only  FS  has  R  >  95% 
with  R  =  86.2%  for  GCI,  R  =  92.1%  for  GCIc,  91.5%  for  CF,  and  95.7%  for  FS  methods. 
For  0  <  CF  <  0.9,  all  methods  show  R  >  93%.  However  for  18%  of  the  data  (CF  =  1.1- 
2.0),  only  47.4%,  71.2%,  72.9%  confidence  interval  can  be  achieved  for  GCI,  GCIc,  and 

FS,  respectively,  while  R  =  89.3%  for  the  FS  method.  For  the  range  of  CF  (0-2),  all  the 
methods  meet  X-tS-  >1.2 .  However,  when  the  CF  (0-2)  is  broken  into  different  ranges, 

only  FS  method  always  satisfies  X-tSx.  >1.2.  For  GCI,  GCIc,  and  CF  methods,  CF 
ranges  where  X-tS-<  1.2  are  (1. 1-2.0)  for  GCI,  (1. 1-1.6)  for  GCIc,  and  (0.9-1. 6)  for  the 

CF  method.  This  puts  a  very  high  risk  of  using  the  three  methods  on  practical 
applications  since  many  verification  studies  for  complex  geometries  and  high  Reynolds 
number  flows  often  approach  the  AR  oscillatorially  with  the  different  sets  of  grids,  which 
has  been  demonstrated  in  many  previous  simulations  in  CFD  workshops  and  an  example 
for  ship  hydrodynamics  as  provided  in  the  next  section.  It  should  be  also  noted  that 
different  choices  of  ranges  of  CF  values,  for  example  0-0.4,  0-0.8,  0.8- 1.2,  1.2- 1.6,  and 
1. 6-2.0,  only  cause  minor  changes  of  the  fraction  bounded  and  the  conclusions  drawn 
above  are  general. 

Table  4  additionally  provides  statistics  for  the  16  CF  where  at  least  4  FSa  are 
available  after  removing  the  outliers.  For  the  1 1  CF  that  are  less  than  1,  all  methods  have 
R  =  100%.  Near  the  AR  where  CF  =  1.005,  GCI,  GCIc,  and  FS  have  the  same  reliability 
R  =  96%,  which  is  larger  than  R  =  94%  for  the  CF  method.  For  CF  close  to  the  AR  but 
larger  than  1  (CF  =  1.015  and  1.095),  R=  100%  for  all  methods.  However,  at  CF  =  1.145, 
only  R  =  69.23%  is  achieved  for  GCI  and  CF  methods,  R  =  76.92%  for  GCIc  method,  but 
with  92.31%  for  the  FS  method.  When  solutions  are  far  away  from  the  asymptotic  range 
(CF  =  1.295),  GCI  only  provides  R  =  75%  while  the  other  three  methods  can  achieve  R  = 
100%.  When  CF  <  1,  the  mean  FSa  X  is  the  same  for  the  GCI  and  GCIc  methods.  For 
CF  =  0.625,  0.675,  and  0.745,  X(GCI)  =  X(GCIc)  <  X(FS)  <  X(CF) .  For  CF  =  0.825 

and  0.865,  X(GCI)  =  X(GCIc)  <  X(CF)  <  X(FS) .  For  CF  =  0.895,  0.915,  0.945,  0.955, 
0.965,  and  0.985,  X(CF)<X(GCI)  =  X(GCIc)<X(FS).  For  CF  =  1.005,  1.015,  and 
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1.095,  ~X(CF)<Jc(GCI)<X(GGIc)<Jc(FS).  For  CF  =  1.145  and  1.295, 

X(GCI)<X(CF)<X(GCIc)<X(FS).  The  different  magnitudes  of  X  follow  the 
theoretical  values  of  FS  at  different  CF  as  shown  in  Figure  1.  Since  a  larger  FS  will  result 
in  a  larger  standard  deviation,  the  s_  follow  the  correlations  above  for  X  at  the  same  CF. 

FS  is  the  only  method  that  satisfies  X-tS->  1.2  at  all  the  16  CF.  X-tS->  1.2  is 

observed  for  GCI  at  CF  =  0.915,  1.145,  and  1.295;  for  GCIc  at  CF  =  0.915,  and  1.145; 
and  for  CF  at  CF  =  0.915,  0.945,  0.955,  0.965,  0.985,  1.005,  1.015  and  1.145.  This 
indicates  that  the  GCI  and  GCIC  don’t  have  sufficient  FS  far  away  from  the  asymptotic 
range  for  CF  >  1,  whereas  CF  method  does  not  provide  enough  FS  close  to  the  AR. 

VII.  EXAMPLE  FOR  SHIP  HYDRODYNAMICS  APPLICATIONS 

To  evaluate  the  behaviors  of  the  four  methods  for  practical  applications,  they  are 
applied  for  a  recent  study  [32]  that  used  computational  towing  tank  procedures  for  single 
run  curves  of  resistance  and  propulsion  for  the  high-speed  transom  ship  Athena  barehull 
with  a  skeg  using  the  general-purpose  solver  CFDShip-Iowa-V.4  [33].  In  this  study,  Xing 
et  al.  investigated  the  issue  of  achieving  the  AR  by  continuously  refining  the  grid  from 
the  coarsest  grid  (grid  7  with  360,528  points)  to  the  finest  grid  (grid  1  with  8.1  million 
points)  for  the  Athena  bare  hull  with  skeg  with  2  degrees  of  freedom  (pitch  and  heave)  at 
Froude  number  (Fr)  0.48.  The  grids  are  designed  with  a  systematic  grid  refinement  ratio 
r  =  2" 25 ,  which  allows  9  sets  of  grids  for  verification  and  validation  (V&V)  with  5  sets 
with  r  =  2025  (5,  6,  7;  4,  5,  6;  3,  4,  5;  2,  3,  4;  and  1,  2,  3),  3  sets  with  r  =  205(3,  5,  7;  2,  4, 
6;  and  1,  3,  5),  and  1  set  with  r  =  2° 75  (1,4,7).  The  distribution  of  iterative  errors 
0.15^,  <  U ,  <0.3 Sfme  for  grids  1  to  7  is  shown  in  Figure  5(b)  and  5(d)  for  resistance 

coefficients  and  motions,  respectively.  U}  is  of  the  same  order  of  magnitude  for  all  the 
grids,  which  suggests  that  it  is  mainly  detennined  by  the  iterative  method  applied  and 
independent  of  grid  resolutions.  As  shown  in  Table  5,  CTX  monotonically  converges  for 

grids  (2,  4,  6),  (1,  3,  5),  (4,  5,  6),  (3,  4,  5),  (2,  3,  4),  and  (1,  2,  3),  of  which  grids  (1,  2,  3) 
have  the  smallest  grid  uncertainty  and  grids  (3,  4,  5)  are  closest  to  the  AR  based  on  CF 
closest  to  one.  CTX  oscillatorially  diverges  on  grids  (5,  6,  7),  and  monotonically  diverges 
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on  grids  (3,  5,  7)  and  (1,  4,  7).  All  the  diverged  solutions  involve  the  coarsest  grid  7, 
which  is  likely  due  to  the  insufficient  resolution  of  the  coarsest  grid  7.  As  shown  in 
Figure  5(a),  CTX  for  grid  7  does  not  follow  the  trend  as  shown  for  grids  6-1.  Figure  5(a) 
also  shows  frictional  and  pressure  resistance  coefficients  Cfx  and  Cpx  on  all  the  grids. 
Figure  5(b)  shows  the  magnitudes  of  the  relative  changes  of  solutions  eN  between  two 
successive  grids  with  respect  to  the  solutions  on  the  finest  grid  1 .  When  grids  are  refined 
from  5  to  1,  eN  systematically  decreases  for  CTX  and  C/Y  while  oscillatory  decreases  for 

Cpx.  UG  for  grids  (4,  5,  6)  is  unreasonable  large  as  it  is  too  far  away  from  the  AR.  Ug 

for  grid  studies  (2,  3,  4)  and  (1,  2,  3)  using  GCI,  GClc,  and  CF  are  unreasonable  small 
due  to 

Table  5.  Verification  study  for  Ctx  of  Athena  bare  hull  with  skeg  (F/^0.48)*. 


Grids 

r 

*0 

Po 

CF 

UG  (%) 

GCI 

GClc 

CF 

FS 

ra 

0.5 

2 

0.63 

1.32 

0.58 

3.34 

3.34 

4.90 

4.17 

1,3,5 

0.5 

2 

0.40 

2.66 

1.51 

0.72 

1.09 

1.16 

2.90 

4,5,6 

0.25 

2 

0.97 

0.16 

0.07 

52.7 

52.7 

125.2 

85.2 

ESS 

0.25 

2 

0.80 

1.27 

0.59 

4.98 

4.98 

7.23 

6.20 

2,3,4 

0.25 

2 

0.60 

2.98 

1.64 

1.07 

1.75 

1.95 

6.64 

1,2,3 

0.25 

2 

0.50 

4.00 

2.42 

0.58 

1.40 

1.11 

- 

*Ua  is  %Sfme;  Cn-is  based  on  static  wetted  area; 


Table  6.  Verification  study  for  motions  of  Athena  bare  hull  with  skeg  (AG^O.48)* 


*UG  is  %Sfine 
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Figure  5.  Verification  for  resistance  and  motions  for  Athena  bare  hull  with  skeg 

(F>^0.48):  (a)  resistance  coefficients,  (b)  relative  change  eN  =  |(.S'V ,  -  .S',  )/. S',  xlOO 

and  iterative  errors  for  resistance  coefficients,  (c)  sinkage  and  trim,  (d)  relative 
change  sn  and  iterative  errors  for  sinkage  and  trim. 


the  deficiency  discussed  before  while  Ug  of  the  FS  method  for  these  two  studies  are  more 
reasonable  on  (2,  3,  4)  but  unfortunately  invalid  for  (1,  2,  3)  due  to  CF  >  2  caused  by  the 
contamination  of  the  iterative  error  on  the  fine  grid.  As  shown  by  Figure  5(b),  separating 
iterative  errors  from  grid  uncertainties  is  problematic  for  the  finer  grids  since  iterative 
and  grid  uncertainties  are  of  the  same  order  of  magnitude.  Implementation  of  more 
accurate  and  efficient  iterative  methods  to  speed  up  the  convergence  (e.g.,  multigrid)  will 
be  necessary.  Flowever,  Sn  of  the  current  study  does  show  systematic  decreasing  for  Ctx 
and  Cfx  and  oscillatory  decreasing  for  CPx.  Ctx,  Cfx,  and  Cpx  show  different  rates  of 


24 


approaching  the  AR.  Table  6  shows  that  motions  are  more  difficult  to  converge  than 
resistance  coefficients.  Mono  tonic  converged  solutions  are  achieved  only  for  trim  on 
grids  (1,  3,  5)  and  (4,  5,  6).  For  trim  on  grids  (1,  3,  5),  FS  method  again  predicts  more 
reasonable  Ug  than  the  other  three  methods.  For  trim  on  grids  (4,  5,  6),  FS  method  also 
predicts  larger  Ug  than  GC1  and  GC1C  methods  but  lower  Ug  than  the  CF  method  that 
seems  to  be  too  conservative  at  CF  =  0.4. 

VIII.  CONCLUSIONS 

A  factor  of  safety  method  for  quantitative  estimates  of  grid  and  time-step 
convergence  uncertainties  for  CFD  solutions  is  derived  to  remove  the  two  deficiencies  of 
using  GCI,  GCIc,  and  CF  methods,  i.e.,  unreasonably  small  uncertainty  when  CF  >  1  and 
lack  of  statistical  analysis  to  prove  95%  confidence  for  the  estimated  uncertainties  to 
bound  the  true  error.  The  approach  follows  the  CF  method  but  by  reflecting  uncertainty 
instead  of  FS  for  CF  <  1  for  CF  >  1  with  respect  to  the  AR  (CF  =1).  Additionally,  the  FS 
method  provides  flexibility  for  factors  of  safety  at  CF  =  0  and  CF  =  1,  which  are 
determined  by  the  overall  reliability  of  the  method  and  lower  band  of  the  confidence 
interval  of  the  true  mean  based  on  statistical  analysis.  The  statistical  analysis  is  based  on 
a  large  sample  of  AB  and  NB  covering  17  studies  with  96  variables  and  304  individual 
grid  triplet  studies,  as  well  as  one  practical  application  in  ship  hydrodynamics.  Results 
showed  that  FS  method  is  the  only  one  that  provides  at  least  95%  confidence  that  the 
uncertainty  estimate  will  bound  the  true  error  (FSa  >  1)  for  the  304  grid  convergence 
studies,  i.e.,  confidence  intervals  are  86.2%,  92.1%,  91.5%,  and  95.7%  for  the  GCI, 
GCIc,  CF,  and  the  FS  methods,  respectively.  For  20%  of  the  selected  data  when  1.1  <  CF 
<  2.0,  GCI,  GCIc,  and  CF  fail  as  only  47.4%,  71.2%,  72.9%  confidence  interval  can  be 
achieved  for  them,  respectively,  while  89.3%  is  achieved  for  the  FS  method.  The  FS 
method  is  also  the  only  one  that  has  at  least  95%  confidence  the  lower  band  of  the 
confidence  interval  for  the  true  mean  FSa  is  larger  than  1 .2  for  different  studies,  different 
variables,  different  ranges  of  CF,  and  different  single  CF  values  where  multiple  FSa  are 
available. 

Although  a  large  sample  of  populations  with  either  analytical  or  numerical 
benchmark  data  have  been  used,  additional  rigorous  verification  studies  (AB/NB  or 
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industrial  applications)  are  needed  to  further  validate  the  recommended  FSo  and  FSi, 
especially  those  demonstrating  that  the  AR  is  achieved  when  grids  are  refined.  This  has 

two  benefits:  (1)  reduced  Sx  and  thus  k,  which  will  likely  provide  a  “converged  curve” 
for  the  factor  of  safety  for  different  CF;  (2)  enabled  Chi-square  analysis  for  evaluation  of 
confidence  interval  for  S-.  Further  evaluation  of  the  LSM  is  needed  especially  using 

statistical  analysis  tools  in  this  study.  There  are  other  unresolved  issues  and  complex 
factors  such  as  mixed  numerical  schemes,  coupled  numerical  and  modeling  errors  for 
large  eddy  simulation  and  detached  eddy  simulation  and  single-grid  estimator.  Even 
though  only  a  small  fraction  (5.6%)  of  data  is  outliers,  it  is  worthy  investigating  why. 
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Appendix  1.  Analytical  and  Numerical  Benchmarks 


Study  1  One-Dimensional  Wave  Equation 

As  shown  in  Figure  Al,  with  the  refinements  of  the  grids,  the  errors  monotonically 
decrease  and  solutions  approach  the  AR  as  CF  changes  from  CF  =  0.4  on  the  coarsest 
grid  to  CF  =  1  on  the  finest  grid.  FSa  for  GCI  and  GCIc  are  the  same  as  CF  <  1  and 
smaller  than  that  predicted  by  the  FS  method  except  at  CF  =  1  where  FSa  using  the  three 
methods  are  the  same.  FSa  for  GCI  and  GCIc  are  smaller  and  larger  than  those  predicted 
by  the  CF  method  for  CF  <  0.875  and  0.875  <  CF  <  1,  respectively.  FSa  of  the  FS 
method  are  smaller  and  larger  than  those  predicted  by  the  CF  method  for  CF  <0.8  and 
0.8  <  CF  <  1,  respectively.  These  observations  are  consistent  with  the  FS  shown  in  Figure 
1.  The  average  FSa  are  1.44,  1.44,  1.61,  and  1.60  for  GCI,  GCIc,  CF,  and  FS  methods, 
respectively.  The  standard  deviations  of  FSa  are  0.27,  0.27,  0.78,  and  0.49  for  GCI, 
GCIc,  CF,  and  FS  methods,  respectively. 

Study  2  Two-Dimensional  Laplace  Equation 

As  shown  in  Figure  A2,  with  the  refinements  of  the  grids,  the  errors  monotonically 
decrease.  Solutions  approach  the  AR  (CF  =  1)  monotonically  for  constant  Dirichlet 
boundary  conditions  and  oscillatorially  for  non-constant  value  Dirichlet  boundary 
conditions.  Comparison  of  FSa  magnitudes  for  different  methods  shows  similar  trend  as 
discussed  in  study  1  and  shown  in  Figure  1.  The  average  FSa  are  1.30,  1.31,  1.23,  and 
1.38  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively.  The  standard  deviations  of  FSA 
are  0.04,  0.04,  0.09,  and  0.07  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  3  Two-Dimensional  Driven  Cavity  (Re=l,000) 

As  shown  in  Figure  A3,  only  two  sets  of  grid  study  are  available  and  thus  it  is  difficult  to 
draw  any  conclusions  on  the  convergence  characteristics  when  the  grids  are  refined. 
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However,  CF  is  closer  to  1  and  errors  decrease  with  the  refinements  of  the  grids.  One 
exception  is  for  the  minimum  stream  function  using  the  BS  scheme,  which  shows  a  larger 
Ug  on  the  finer  grid.  Comparison  of  FSa  magnitudes  for  different  methods  shows  similar 
trend  as  shown  in  Figure  1.  All  verification  methods  except  the  FS  method  predict  FSa  < 
1  for  CF  >  1.  The  average  FSa  are  1.12,  1.16,  1.23,  and  1.58  for  GCI,  GCIc,  CF,  and  FS 
methods,  respectively.  The  standard  deviations  of  FSa  are  0.41,  0.37,  0.51,  and  0.62  for 
GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  4  Two-Dimensional  Natural  Convection  Flows  in  Square  Cavities  (Ra=104) 

As  shown  in  Figure  A4,  all  variables  have  3  grid-triplet  studies  except  Umon  that  has  only 
two  grid-triplet  study.  When  grids  are  refined,  all  errors  monotonically  decrease  and  all 
variable  solutions  either  linearly  or  oscillatorially  approach  the  AR  (CF  =1).  Comparison 
of  FSa  magnitudes  for  different  methods  shows  similar  trend  as  shown  in  Figure  1.  FSa 
oscillatorially  decrease  when  CF  is  approaching  the  AR  from  CF  <  1.  For  CF  >  1,  GCIc 
predicted  higher  FSa  than  GCI  and  CF  methods.  For  CF  >  1.16,  uncertainty  estimates 
using  GCI  do  not  bound  the  true  error  (FSa  <  1)  and  GCIc  and  CF  have  similar 
magnitudes  of  FSa  around  1 .4  while  the  FS  method  has  much  larger  FSa.  The  average 
FSa  are  1.49,  1.55,  1.75,  and  1.90  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively.  The 
standard  deviations  of  FSa  are  0.26,  0.25,  0.68,  and  0.45  for  GCI,  GCIc,  CF,  and  FS 
methods,  respectively. 

Study  5  Two-Dimensional  Natural  Convection  Flows  in  Square  Cavities  (Ra=105) 

As  shown  in  Figure  A5,  all  variables  have  4  grid-triplet  studies.  When  the  grids  are 
refined,  all  errors  monotonically  decrease  except  the  two  coarsest  sets  of  grids  for  Umax 
and  all  variable  solutions  either  monotonically  or  oscillatorially  approach  the  AR.  Similar 
to  the  previous  studies,  uncertainty  estimates  using  GCI,  GCIc,  and  CF  methods  do  not 
bound  the  true  error  for  some  points  at  CF  >  1  (FSa  <  1)  while  the  FS  method  always 
bounds  the  true  error  (FSa  <  1).  The  average  FSA  are  1.84,  1.96,  2.64,  and  2.74  for  GCI, 
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GCIc,  CF,  and  FS  methods,  respectively.  The  standard  deviations  of  FSa  are  0.74,  0.73, 
1.77,  and  1.48  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  6  Two-Dimensional  Natural  Convection  Flows  in  Square  Cavities  (Ra=106) 

As  shown  in  Figure  A6,  the  number  of  grid-triplet  studies  is  5  for  all  variables  except  4 
for  Umax  and  Numax  and  2  for  Vmax.  When  the  grids  are  refined,  all  errors  monotonically 
decrease  and  all  variable  solutions  either  monotonically  or  oscillatorially  approach  the 
AR.  Similar  to  the  previous  studies,  GCI,  GCIc,  and  CF  methods  do  not  bound  the  true 
error  for  some  points  at  CF  >  1  while  the  FS  method  always  bounds  the  true  error.  The 
average  FSa  are  1.47,  1.53,  1.91,  and  1.89  for  GCI,  GCIc,  CF,  and  FS  methods, 
respectively.  The  standard  deviations  of  FSA  are  0.40,  0.40,  0.99,  and  0.64  for  GCI, 
GCIc,  CF,  and  FS  methods,  respectively. 

Study  7  Two-Dimensional  Backward-Facing  Step 

As  shown  in  Figure  A7,  with  the  refinement  of  the  grids,  the  errors  oscillatorially 
decrease  with  several  orders  of  magnitude  larger  than  those  in  other  studies  and  solutions 
approach  the  AR  with  huge  oscillations  of  CF.  This  suggests  that  the  grids  are  still  far 
from  the  AR  for  such  a  high  Reynolds  number  flows.  When  CF  >  1.4,  GCI  does  not 
bound  the  true  error  while  all  other  methods  bounded  the  true  error  with  FSa  predicted  by 
the  FS  method  is  the  largest.  The  average  FSa  are  0.92,  1.39,  1.49,  and  15.08  for  GCI, 
GCIc,  CF,  and  FS  methods,  respectively.  The  standard  deviations  of  FSa  are  0.33,  0.15, 
0.22,  and  19.12  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  8  Square  Cavity  with  Moving  Top  Wall  (Re=100) 

As  shown  in  Figure  A8,  4  variables  only  have  two  grid-triplet  studies  and  3  variables 
have  3  grid-triplet  studies,  which  suggests  that  further  refinement  of  the  grids  are 
necessary.  However,  all  variable  errors  linearly  decrease  when  the  grids  are  refined. 
When  CF  >  1.4,  all  methods  except  GCI  predict  enough  uncertainty  estimates  to  bound 
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the  true  errors.  The  average  FSa  are  1.56,  1.96,  2.20,  and  6.27  for  GCI,  GCIc,  CF,  and  FS 
methods,  respectively.  The  standard  deviations  of  FSa  are  0.48,  0.45,  0.78,  and  6.38  for 
GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  9  Square  Cavity  with  Moving  Top  Wall  (Re=l,000) 

As  shown  in  Figure  A9,  all  variables  only  show  monotonic  convergence  on  two  grid- 
triplet  studies  except  3  grid  studies  for  the  x  velocity  using  UDS  at  0  degree.  When  the 
grids  are  refined,  errors  for  UDS-x  velocity-0  degree  oscillatorially  decrease  while  errors 
for  all  other  variables  linearly  decrease.  When  CF  =  1  and  1.15,  all  methods  fail  to  bound 
the  true  error  except  the  FS  method.  The  average  FSa  are  2.42,  2.46,  4.12,  and  3.32  for 
GCI,  GCIc,  CF,  and  FS  methods,  respectively.  The  standard  deviations  of  FSA  are  1.55, 
1.52,  3.72,  and  2.41  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  10  Cubic  Cavity  with  Moving  Top  Wall  (Re=100) 

As  shown  in  Figure  A10,  only  two  grid-triplet  studies  are  available.  However,  all  errors 
linearly  decrease  when  the  grids  are  refined.  Uncertainty  estimates  predicted  by  all 
methods  bound  the  true  errors.  The  average  FSa  are  2.76,  2.90,  3.93,  and  4.10  for  GCI, 
GCIc,  CF,  and  FS  methods,  respectively.  The  standard  deviations  of  FSa  are  1.43,  1.28, 
2.72,  and  1.61  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  11  Axisymmetric  Turbulent  Flow  through  a  Valve 

As  shown  in  Figure  All,  similar  to  study  10,  only  two  grid-triplet  studies  are  available 
and  thus  no  conclusive  behavior  of  the  variables  to  approach  the  AR  can  be  drawn. 
However,  all  errors  linearly  decrease  when  the  grids  are  refined.  Uncertainty  estimates 
predicted  by  all  methods  bound  the  true  errors.  The  average  FSa  are  3.22,  3.22,  5.89,  and 
4.46  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively.  The  standard  deviations  of  FSA 
are  2.58,  2.58,  5.14,  and  3.73  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 
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Study  12  One-Dimensional  Steady-State  Convection-Diffusion  Process  without 
Source  Term 

As  shown  in  Figure  A 12,  all  variables  have  four  grid-triplet  studies.  When  the  grids  are 
refined,  errors  for  all  variables  decrease  and  solutions  are  approaching  the  AR. 
Uncertainty  estimates  of  all  methods  bound  the  true  errors  except  at  CF  =  0.91  where 
only  the  FS  method  bound  the  true  error.  The  average  FSa  are  1.50,  1.51,  1.51,  and  1.63 
for  GCI,  GCIc,  CF,  and  FS  methods,  respectively.  The  standard  deviations  of  FSa  are 
0.35,  0.34,  0.41,  and  0.38  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  13  Heat  Transfer  from  Isothermal  Cylinder  Enclosed  by  a  Square  Duct 

As  shown  in  Figure  A 13,  only  two  grid-triplet  studies  are  available  for  all  the  variables 
except  3  grid-triplet  studies  for  the  UDS -temperature.  When  the  grids  are  refined,  all 
errors  linearly  decrease  except  the  uncertainty  estimates  using  GCI,  GCIc,  and  CF 
methods  for  temperature  with  the  SMART  scheme  increase.  When  1.2  <  CF  <  1.4, 
uncertainty  estimates  predicted  by  all  the  methods  do  not  bound  the  true  errors  except  the 
FS  method.  The  average  FSa  are  1.39,  1.63,  1.65,  and  2.50  for  GCI,  GCIc,  CF,  and  FS 
methods,  respectively.  The  standard  deviations  of  FSa  are  0.47,  0.47,  0.33,  and  0.81  for 
GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  14  Premixed  Methane/ Air  Laminar  Flat  Flame  on  a  Perforated  Burner 

As  shown  in  Figure  A14,  there  are  5  and  4  grid-triplet  studies  for  UDS  and  SMART 
methods,  respectively.  When  the  grids  are  refined,  errors  for  most  variables  continuously 
decrease,  solutions  using  the  UDS  method  approach  the  AR  either  linearly  or 
oscillatorially,  and  solutions  using  the  SMART  method  show  oscillations  far  away  from 
the  AR.  For  CF  close  to  1,  GCI,  GCIc  and  CF  predict  FSa  around  1.  For  CF  >  1,  FSa 
predicted  by  the  three  methods  show  large  oscillations  of  which  some  values  are  below  1 . 
The  FS  method  always  provides  sufficient  large  FSa-  However,  there  is  one  grid-triplet 
study  at  CF  =  0.51  where  none  of  the  methods  can  bound  it.  The  reason  behind  this  is 
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unknown.  The  average  FSa  are  1.63,  1.66,  2.13,  and  2.25  for  GCI,  GCIc,  CF,  and  FS 
methods,  respectively.  The  standard  deviations  of  FSa  are  0.35,  0.35,  0.87,  and  1.02  for 
GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  15  Data  for  “Exact”  Grid  Convergence  Set 

As  shown  in  Figure  A15,  solutions  for  the  “exact”  grid  convergence  study  are  always  in 
the  AR  when  there  is  no  perturbation  or  with  perturbation  #4.  Other  cases  only  have  two 
grid-triplet  studies.  Overall  errors  decrease  when  the  grids  are  refined.  For  CF  >  1  when 
solutions  are  far  away  from  the  AR,  FSa  for  the  FS  method  increases  while  FSa  for  other 
methods  decrease.  The  average  FSa  are  1.79,  1.81,  2.24,  and  2.18  for  GCI,  GCIc,  CF,  and 
FS  methods,  respectively.  The  standard  deviations  of  FSA  are  0.85,  0.88,  1.76,  and  1.44 
for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  16  Beam  Bending  Problem  for  2nd  Series  on  Grid  Convergence 

As  shown  in  Figure  A 16,  there  are  3  grid-triplet  studies  for  the  three  codes  investigated. 
When  the  grids  are  refined,  both  the  errors  and  the  correction  factors  show  oscillations. 
Uncertainty  estimates  predicted  by  all  methods  fail  to  bound  the  true  errors  when  0.6  < 
CF  <  1.2.  The  average  FSa  are  0.59,  0.64,  1.00,  and  0.96  for  GCI,  GCIc,  CF,  and  FS 
methods,  respectively.  The  standard  deviations  of  FSa  are  0.068,  0.061,  0.228,  and  0.095 
for  GCI,  GCIc,  CF,  and  FS  methods,  respectively. 

Study  17  Beam  Bending  Problem  for  3rd  Series  on  Grid  Convergence 

As  shown  in  Figure  A17,  there  are  only  two  grid-triplet  studies  for  the  three  codes 
investigated.  When  the  grids  are  refined,  the  errors  decrease.  Uncertainty  estimates 
predicted  by  all  methods  fail  to  bound  the  true  errors  when  CF  >  1.3.  The  average  FSa 
are  0.97,  1.03,  1.69,  and  1.52  for  GCI,  GCIc,  CF,  and  FS  methods,  respectively.  The 
standard  deviations  of  FSa  are  0.58,  0.52,  1.31,  and  0.68  for  GCI,  GCIc,  CF,  and  FS 
methods,  respectively. 
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Figure  Al: 


(a)  (b) 

ID  Wave  Equation  [6]:  (a)  error,  uncertainties,  and  correction  factor, 
(b)  actual  factor  of  safety 
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Figure  A2:  2D  Laplace  Equation  [5,7]:  (a)  non-constant  value  Dirichlet  boundary 
conditions  (Eca  and  Hoekstra,  2000),  (b)  non-constant  value  Dirichlet 
boundary  conditions  (Iowa  recalculated),  (c)  constant  value  Dirichlet 
boundary  conditions  (Iowa  calculation),  (d)  actual  factor  of  safety. 
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Figure  A3:  2D  Driven  Cavity  [20,  21]:  (a)  maximum  stream  function  using  BS  scheme, 
(b)  maximum  stream  function  using  upwind  3  scheme,  (c)  maximum  stream 
function  Kawamura  scheme,  (d)  maximum  vorticity  using  BS  scheme, 

(e)  maximum  vorticity  using  upwind-3  scheme,  (f)  maximum  vorticity  using 
Kawamura  scheme,  (g)  minimum  stream  function  using  BS  scheme,  (h)  actual 
factor  of  safety. 
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Figure  A4:  2D  Natural  Convection  Flows  in  Square  Cavities  at  Ra=  1 04  [22]:  (a)  Umon, 
(b)  Vmon,  (c)  Tmon,  (d)  Umax,  (e)  Vmax,  (f)  Numax,  (g)  Nu,  (h)  actual  factor  of 
safety. 
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Figure  A6:  2D  Natural  Convection  Flows  in  Square  Cavities  at  Ra=  1 0fi  [22]:  (a)  Umon, 
(b)  Vmon,  (c)  Tmon,  (d)  Umax,  (e)  Vmax,  (f)  Numax,  (g)  Nu,  (h)  actual  factor  of 
safety. 
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Figure  A7:  2D  Backward-facing  Step  at  Re=1.5xl05  [23,  24,  25]:  (a)  error,  uncertainties, 
and  correction  factor,  (b)  actual  factor  of  safety 
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Figure  A8:  Square  cavity  with  moving  top  wall,  Re=100  [26,  27]:  (a)  UDS-x  velocity-0 
deg;  (b)  UDS-x  velocity-60  deg;  (c)  UDS-y  velocity-0  deg;  (d)  SMART-x 
velocity-0  deg;  (e)  SMART-x  velocity-60  deg;  (f)  SMART-y  velocity-0  deg; 
(g)  SMART-y  velocity-60  deg;  (h)  actual  factor  of  safety. 
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Figure  A9:  Square  cavity  with  moving  top  wall  at  Re=1000  [26,  27]:  (a)  UDS-x  velocity- 
0  deg;  (b)  UDS-x  velocity-60  deg;  (c)  UDS-y  velocity-0  deg;  (d)  UDS-y 
velocity-60  deg;  (e)  SMART-x  velocity-0  deg;  (f)  SMART-x  velocity-60  deg; 
(g)  SMART -y  velocity-0  deg;  (h)  SMART-y  velocity-60  deg;  (i)  actual  factor 
of  safety. 
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Figure  A 10:  Cubic  cavity  with  moving  top  wall,  Re=100  [26,  27]:  (a)  x-velocity-UDS 
(b)  y-velocity-UDS;  (c)  z-velocity-UDS,  (d)  actual  factor  of  safety. 


Grid  refinement  ratio  (Ax/Ax^) 


1  1.2  1.4  1.6  1.8  2 

Grid  refinement  ratio  (Ax/Ax^J 


CF 


(e) 

Figure  All:  Axisymmetric  turbulent  flow  through  a  valve  [28,  29]:  (a)  TKE-PLDS; 

(b)  dissipation  rate  of  TKE-PLDS;  (c)  radial-velocity-SMART;  (d)  axial- 
velocity-SMART;  (e)  actual  factor  of  safety. 
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Figure  A 12:  One-dimensional  steady-state  convection-diffusion  process  without  source 
term,  with  constant  transport  properties  and  with  Dirichlet  boundary 
conditions,  Pe=l  and  Pe=10  [1 1]:  (a)  UDS-Pe=l,  (b)  UDS-Pe=10, 

(c)  SMART-Pe=l,  (d)  SMART-Pe=10,  (e)  actual  factor  of  safety. 
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Figure  A13:  Heat  transfer  from  an  isothermal  cylinder  enclosed  by  a  square  duct  [27,  28, 
29]:  (a)  UDS-y-velocity,  (b)  UDS-temperature,  (c)  SMART-x-velocity, 

(d)  SMART -y-velocity,  (e)  SMART-temperature,  (f)  actual  factor  of  safety. 
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Figure  A14:  Premixed  methane/air  laminar  flat  flame  on  a  perforated  burner  [28,  30,  31] 
(a)  UDS-radial-velocity,  (b)  UDS-axial-velocity,  (c)  UDS-temperature, 

(d)  UDS-CH4,  (e)  UDS-02,  (f)  UDS-C02,  (g)  UDS-H20,  (h)  SMART- 
radial-velocity,  (i)  SMART-axial-velocity,  (j)  SMART-temperature, 

(k)  SMART-CH4,  (1)  SMART-02,  (m)  SMART-C02,  (n)  SMART  H20, 
(o)  actual  factor  of  safety. 
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Figure  A15:  Data  for  “exact”  grid  convergence  set  [8]:  (a)  exact  to  F=600,  (b)  perturbed 
#3,  (c)  perturbed  #4,  (d)  all±5,  (e)  actual  factor  of  safety. 
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Figure  A 16:  Beam  bending  problem  for  2nd  series  on  grid  convergence  [8]:  (a)  beam 
bending  stress  (Code  1),  (b)  beam  bending  stress  (Code  X),  (c)  beam 
bending  stress  (Code  W),  (d)  actual  factor  of  safety. 
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Figure  A 17:  Beam  bending  problem  for  3rd  series  on  grid  convergence  [8]:  (a)  beam 
bending  stress  (Code  1,  code  X),  (b)  beam  bending  stress  (Code  W), 

(c)  Beam  end  deflection  (Code  1),  (d)  Beam  end  deflection  (Code  X), 
(e)  actual  factor  of  safety. 
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